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ABSTRACT 


A numerical filtering technique useful in removing the diurnal 
component from surface data of magnetic field measurements is de- 
scribed. Derivations of formulas used to compute filter weights and 
several examples of the current use of such filters are presented and 
outlined. 



CONTENTS 

Abstract ii 

INTRODUCTION 1 

THE LEAST SQUARES APPROXIMATION TO T(f) FOR 

A LOW PASS FILTER 8 

THE CHEBYSHEV APPROXIMATION TO T(f) 

FOR A LOW PASS FILTER 18 

SCALING FILTERS AND SHIFTING FOR 

BANDPASS RESPONSE 23 

CONCLUDING REMARKS 33 

References 35 


THE DESIGN OF NUMERICAL FILTERS 
FOR GEOMAGNETIC DATA ANALYSIS 


by 

Kenneth W. Behannon and Norman F. Ness 
Goddard Space Flight Center 


INTRODUCTION 

In investigating the transient time variations that occur throughout the geomagnetosphere, 
effective comparison can be made between the data representing the magnetic field measurements 
from surface observatories and satellite magnetometer data only if the periodic effects due to the 
earth r s rotation, such as Sq , are removed from the ground observatory data. One numerical 
technique which has proven useful in removing the diurnal component from surface data is that of 
numerical filtering. It is the purpose of this discussion to describe briefly this process and to 
show how it has been developed and successfully applied. Although the applications described 
here are rather limited in scope, techniques such as these lend themselves so well to automatic 
data processing that they find general application in many studies of geophysical time or space 
series (Reference 1). 

For processing data in such a way as to selectively remove certain periodic components, one 
can construct linear operators that produce the same effect when applied to experimental data as 
that produced by electrical and optical filters. Analogously then, in discussing these numerical 
filters, we shall speak of the input function l(t), the output function O(t), the theoretical transfer 
function T(f) used in the design of the filter, the gain or frequency response W(f) (which is the 
transfer function of the resultant numerical filter), and the phase shift cp(f) of the filter. 

Two important properties of a linear filter are (1) that the output is a linear function of the 
input, i.e., if two inputs Ij(t) and l 2 (t) give the outputs O t (t) and 0 2 (t) , then input 
I 3 (t) = Xj ( t ) + I 2 (t) will give the output 0 3 (t) = O^t) + o 2 (t) ; and (2) its response is independ- 
ent of the time origin, i.e., if an input l(t) results in an output 0(t) , then an input i(t + t 0 ) gives 
an output 0( t + 1 0 ) . 

Most time series of interest in geophysics can be considered to be quasi- stationary time 
series. A stationary time series is a random function of time which may also be a function of 
initial conditions but whose average probability distributions are independent of time (Reference 2). 
The fundamental principles of stationary time series smoothing have been elaborated by Wiener 
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(Reference 3) and others. In application, simple smoothing methods are included among the basic 
techniques of numerical analysis (Reference 4 and Reference 5). Since such data smoothing is 
actually a type of time or space series filtering, it will be instructive to delay further discussion 
of numerical filter design until a brief description of a simple smoothing operator or, equivalently, 
a filter has been given. 

A time series g (t) of equally- spaced data values at intervals At can be smoothed "by threes" 
using the linear formula 


g k _i(t) + g k (t) + g k+1 (t) 


One can see the result of this operation on the frequency spectrum of g(t) by looking at its effect 
on an individual pure sinusoidal time function. 

As an example let 

g k (t) = Real Part [e i " k ^ t ], 
where ^ « 2rrf (frequency) and kAt = t. Then 


-i 


-i)At + e i^kAt + e i<»(k+x)At| 


3 


^ e i<» kA t 


e -i<wAt + e i<u k At 


+ e ia,kAt 



: I e i“kAt (1 + 2 cos <yAt) 


= Si 


/l + 2 cos a)At \ 

l 5 7 


Thus the individual values in the time series are modified by a factor, the transfer function, 
which is independent of k. In the frequency domain, the spectrum will be amplitude modulated by 
the factor (l/3)(l + 2 cos oAt) . Plotting this quantity as a function of frequency reveals that most 
of the frequencies in the upper half of the spectrum are suppressed by this simple smoothing 
process (for a certain limited range of frequencies). In this case the time series has been "low 
pass filtered" by taking weighted running means with all weights identical and equal to 1/N, where 
N is the number of data points used in computing the mean. If "aliasing” of the data occurs (see 
following paragraphs) then contributions from frequencies near f = N/At can occur, where N is an 
integer. 
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In general a numerical filter consists of a set of "weights" w k which determine the actual 
transfer function W(f) of the filter. The design of a numerical filter begins with establishing the 
shape of the data window in the frequency domain which will give the desired effect. Having speci- 
fied the theoretical transfer function, the remainder of the problem consists of determining the 
weights W k in such a way that the actual transfer function, or frequency response, approximates 
the desired one as well as possible. A perfect low pass filter, for example, would leave unaltered 
all frequency components from f = o to the desired cutoff frequency f L and then would suppress 
all frequencies greater than f L . The response of an actual numerical filter can only approximate 
this ideal behavior, with the accuracy of the approximation depending on the values of various 
design parameters. 

As in the simple smoothing process, a numerical filter is applied in such that 


= 



w k y (t + kAt). 


(i) 


The filtering is accomplished by "sliding" the filter along the data, applying it to m + 1 + N data 
points to produce the filtered equivalent of the data point which has been multiplied by W 0 and then 
moving each weight to the next point in the series and repeating the application. Repetition of the 
process until all the data in a given run have been covered produces a series of filtered data points 
which defines the output function 0(t) . Within the precision of the filter these points will trace out 
the input function i(t) with the unwanted high frequency components removed (if a low pass filter 
is being used). 

When experimental data are derived by discretely sampling some phenomenon at equally 
spaced intervals of time, the problem of aliasing may occur in which the sampling rate is low 
enough to confuse two or more frequencies in the data. The net result is that they appear to be the 
same frequency (Figures la and lb). To avoid this problem and hence to define a unique input 
function as described by a set of data points, one must be able to assume that the phenomenon 
studied is spectrally limited to the range |f| < f c , where f c = f s /2, f s being the sampling fre- 
quency and f c being the cut-off or Nyquist* frequency. If such an assumption is valid, then the 
function has been sampled frequently enough so that all significant frequency components are de- 
terminable. This is a result of the sampling theorem of information theory (Reference 2). The 
sampling theorem states that if a function G(t) contains no frequencies higher than w cycles per 
second, then it is completely determined by giving its ordinates at a series of points spaced 1/2W 
seconds apart, the series extending throughout the entire time domain. There is an equivalent 
theorem for the frequency domain. 

We shall now very briefly present a few of the analytical considerations underlying numerical 
filter design. First of all, two mathematical concepts basic to all time series analysis are the 
Fourier transform and the operation of convolution (Reference 6). 

* After H. Nyquist of Bell Telephone Laboratories 
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The Fourier transform F(f) of a function 
g ( t ) is defined as 


F(f) = g(t) e' i27rft dt 


-f 

J -oo 


When g(t) is such that 


li 

T 


1 r 

2T J-t 


|g(t )| 2 dt < 


then 


g(t) 


^-0: 


F(f) e 


i 277- ft 


df . 


F(f) and g(t) are referred to as a trans- 
form pair when the above conditions hold. In 
various phases of the analysis of time series 
the operation of Fourier transforming a func- 
tion is required. However, the function may 
not be known analytically but only at certain 
tabulated values. Thus some approximation to 
the above integration must be made with dis- 
crete data. 



EXAMPLE: AT = 1/3 F c = 1 .5 LET F 0 - 0.5 THEN Fj =2.5 


(b) 


FREQUENCY 


2Fc 



Figure 1— Discrete sampling and the aliasing problem: 
(a) An example of the sampling problem in which two 
different sinusoids (F 0 and F t ) have exactly the same 
sample values at the sample points, (b) A line running 
from zero to infinity along which frequencies are repre- 
sented as points (F g = sampling frequency), (c) The same 
line as in (b)but with the Nyquist frequency, F c , placed 
such that F c = F s /2 = l/(2AT)and F, = 2kF c ± F 0 , where 
AT is the period of the sampling frequency and k is an 
integer multiplier and (d) Aliasing is illustrated by the 
line (b) folded back and forth on itself. Subsequent 
folds occur with the same spacing and all points above 
the same point on the frequency axis appear to be the 
same frequency. After sampling, aliased frequencies 
cannot be separated. 


The convolution of the function g(x) with respect to the function f(x) on the interval [a,b] is 
the function 


h(x) 



f(t) g(x-t) dt = f g(t) f(x-t) dt. 


If the interval is extended to include the entire real line then one obtains the function 


H(x) 



f (t) g(x-t) dt . 
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The latter is defined simply as the convolution of the functions f and g (Reference 7). The 
integrand in the convolution integral will have nonzero values only in the region of overlapping of 
the functions f(t) and g(x - t ) . Thus there are cases where h(x) = H(x) for a and b both finite 
numbers. For example, if [a,b] = [o,x] and the integrand is zero for all values of t < 0, all 
nonzero values of the integrand will lie within [o,x] and hence h(x) = H(x). 

Theoretically, the transfer function of a filter is simply the ratio of the output spectrum to the 
input spectrum, 


T(f) = 


0(0 

1(f)’ 


(2) 


or 


OCO = T(f) * 1(f), 

where 1(f) and 0(f) are the Fourier transforms of the input and output functions i(t) andO(t) . 
Multiplication in the frequency domain of functions possessing Fourier transforms corresponds to 
convolution in the time domain of the transforms, and conversely (Reference 6). Thus 


0(t) = 



I (r) T(t-r) dr, 


or, equivalently, 


(3) 


o(t) = 



T(t) I(t-r) dr. 


The Fourier transform of the function T(t) is the theoretical transfer function of the filter: 


T(f) 


/•CD 
J -CO 


T(t) 


dt , 


(4) 


An arbitrary function can always be expressed as the sum of two component functions of which 
one is odd and the other even. Thus we can write the relation 


T(t) = T E + T 0 . 

From the definition of even and odd functions it follows that 

T(-t) = T E -T 0 . 
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Addition and subtraction of these expressions yield 


' T e =-| [T(t) + T(-t)] , 
T 0 =\ [T(t) -T(-t)]. 


Now we can write Equation 4 as 


T(f) = 



(T e + T 0 ) e iwt 


dt , 


(T e + T 0 ) (cos cot + i sin cjt) dt . 


Because of the effect of integrating even and odd functions between symmetric limits this reduces 
to 


T(f) = 2 I T e cos cot dt + 2i I T 0 sin cot dt. 

Jo Jo 

It is seen immediately that for T(f) to be a real number, we must have T d = 0. This in turn re- 
quires that T(t) = T(-t). 

The theoretical transfer function is approximated by the actual transfer function or frequency 
response of a numerical filter: 

N 

W(f ) = ^ W k (t)e i2,rfkAt = T(f). (5) 

k=-M ' ' 


If M = N then one can write 


W(f) 


= w 0 +2^ (W k e i2,TfkAt + w. k e- 


i2n fkA t 


)• 


(6) 


As in the case of continuous functions, we can express w as the sum of even and odd parts: 


W k = W E + w 0 
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and similarly obtain 


W e= j (W k +W_ k ). 

Wo=7(W k -W_ k ). 

Again as in the ideal case one finds that the necessary and sufficient condition for the transfer 
function W(f) to be a real number for any value of f is that the filter be symmetric (W_ k = W k ). 

If the filter is asymmetric (W„ k = -w k ) , the transfer function is a pure imaginary number. 

The complex number representing the transfer function can also be written 

W(f) = G(f ) e iCp < f > (7) 

where G(f) is the gain of the filter and <p(f) is the phase shift that it produces. G(f) must be an 
even function of f and <p(f ) an odd function of f for real input and output functions. The relation 
indicated by Equation 7 may also be written in the form 

W(f) = G(f ) cos <p(f) + i G(f) sin cp(f). (8) 

It has been shown that if the filter is symmetric, W(f) is a real number. From Equation 8 it is 
seen that this requires that q>(f ) be equal to zero or n . Likewise if the filter is asymmetric and 
hence W(f) is a pure imaginary, we must have that <p(f) = *Hr/2 or - 77 / 2 . In other words, besides 
modifying the amplitude of an input frequency component due to the effect of the gain G(f ), a filter 
with W k = w_ k for all k either has no effect on the phase of the input function or shifts it by 77, 
while a filter with -w k = W_ k for all k produces a phase shift of ±77/2. 

One can thus call the even part of the transfer function the "in phase" portion and the odd 
part the "out of phase" portion due to the effect of each when nonzero. In the most general case 
W k = W E + W 0 with neither part equal to zero, i.e., it contains both in-phase and out-of-phase 
portions, is complex, and phase shifts the input function by an amount 0 < cp < 2n . 

From Equation 6 with W 0 = 0 and w_ k = -w k one obtains the formulas for the asymmetric or 
"sine" filter frequency response 


W(f) = 


k=i 


sin 277 f kAt , 


( 9 ) 


Consider as input the complex sine function l(t) = Ae iat with Fourier transform i(f) = A . Sup- 
pose it is desired that the filter output be the derivative of I(t), i.e., that 0(t) = uul(t). Trans- 
forming O(t) reveals that in the frequency domain we must have 0(f) = icuA. Since 0(f) = T(f) • 1(f), 
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it is hence necessary that T(f) = ico = oje i7T /2. Thus in order to perform differentiation the trans- 
fer function of a filter must be such as to produce a gain of a) - 2ni and a phase shift of W2. 
Because it provides the necessary phase shift, the asymmetric or sine type filter may be used as 
a differentiator. 

For most numerical filtering of geomagnetic time series it is desirable only to attenuate 
certain frequency components without altering the phase. Hence of greatest importance is the 
symmetric or "cosine" filter with frequency response 

N 

W(f) = W 0 + 2^]w k cos 27TfkAt. (10) 

k= 1 

This expression may be used to compute the frequency response characteristics for the designed 
filter once the numerical values of the weights w k are known. In the following sections we shall 
discuss two methods for calculating the values of the weights, given the characteristics of the 
theoretical transfer function T(f) . 


THE LEAST SQUARES APPROXIMATION TO T(f) FOR 
A LOW PASS FILTER 

One approach to the approximation of the theoretical transfer function is through application 
of the least squares technique. In the following development of formulas which can be used to 
calculate filter weights, we shall closely follow the discussion of the subject by Martin (Refer- 
ence 8). Instead of using the frequency f, Martin introduces the normalized frequency 
r - f/f s = f/2f c . He then designates r c to represent the ratio of the desired cutoff frequency to 
the sampling frequency. We prefer to use the parameter p = 2r = f/f c for frequency normali- 
zation and p as the cutoff ratio. 

As stated previously, the problem of filter design consists of determining the M + 1 + N 
weights W k such that the actual transfer function of the filter is defined by Equation 5, or, in 
terms of p, 


N 

W(p).£w k e-P. (11) 

k=-M 

approximates the best, in the least squares sense, the desired transfer function. The transfer 
function for a perfect filter may be written in the form 


T(p) = G(p) e i(p <P>, 


( 12 ) 


<p(p) being the phase shift. 


8 



We shall require that the mean square deviation between T(p) and W(p) over a specified 
interval -p' to +p', given by 


I=7T-rf I T(p) -W(p)| 2 dp, 

2 P J-p' 

be minimized by proper choice of the M + 1 + N weights w k . Thus we can write 

p' N 

I = J_ f | G(p) e i<p ( p > W k e i7T kp | 2 dp, 

P J-p' k=— M 

or, since for z* the conjugate of Z, |z | 2 = Z • Z*, 

_ N 

G(p) e - iCp (P> - 2^ W k e -i 77 kp 


( 13 ) 


r = 2 P 


' I = 

* - n ' 


G(p) e iCp < p > - L w k e 1 kp 


k= -M 


dp. 


(14) 


To minimize the function r , the deviation of that function with respect to each w k must be 
zero. In other words we must have 


dp = 0, (15) 


or 


r _ e i?rkp 

j- p # 

N 

G(p')e- i< P <p) - W n e _i7Tnp 

_ e -i?rkp 

N 

G(p) e ic P <p > - ^ W n e i7Tnp 


n = — M _ 


n =— M _ 


p ^ N /%p — 

f W n [e i7r ^ k - n)p + e- i,T(k - n >p] dp - I G(p) j^e‘ ^ 

J-p' n --M *'”P 


kp-Cp(p)] -i[77k! 


> (P)] ^ 


dp = 0. 


This gives us the relation 


J y W n cos 7 r(k - n)p dp = f G(p) cos [ 77 -kp - <p(p)] dp, 

P # n--M •'-P # 

or, by changing the order of integration and summation, 


(16) 
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I i i i ii min i mi nun 


" r r p ‘ 

> W n cos77(k- n)p dp = I G(p)cos [ 7 zkp - cp (p) ] dp. 

n = -M ''-p' •'-P 


Since tt( k - n)p and rrkp - qp(p) are both odd functions of P , their cosines are even functions of 
p and we may hence write the above integrals as 


L •» [ “ 

n = -M •'0 


S7r(k -n)pdp = G(p)cos [77kp - ^(p) ] dp. 


There is one of these equations for each value of k from -m to n, or, for a symmetrical filter, 
from -N to +N . 

It was previously stated that the phenomenon being studied must be spectrally limited to the 
range |f| < f to avoid aliasing problems. Equivalently, we require | p | < 1 , this leading to 


cos 7T (k - n) pdp 


Hence we are left with 


W = G(p) cos [77-kp - <p(p)] dp, 


so that each W is expressed explicitly. We may write this as 


W - G(p) cos [77kp -cp(p) ] dp + G(p) cos [ 77 kp-qp(p)] dp, 


where p is the cutoff ratio. 

For the ideal low pass filter, G(p) = 1 for 0 < p < P andG(p) = 0 for p > P. Hence 


J o 


>s [77k p - <p(p)] dp, 
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or, for zero phase shift 


■f 

J o 


COS 77 kp dp. 


(23) 


For k = 0 this gives us 


w 


o 



P, 


(24) 


and for k / 0 we have 


\ = 


sin 7 rk P 
77k 


(25) 


Equations 24 and 25 may be used to compute low pass filter weights for sharp cutoff, but they 
lead to an approximation of the ideal transfer function which exhibits a large overshoot for values 
of P slightly smaller or greater than P . This is a manifestation of the Gibbs phenomenon dis- 
cussed in most works on Fourier analysis. This phenomenon occurs near a discontinuity in a 
function which is being approximated by a finite series of size n. As n increases, the position 
at which the maximum occurs moves nearer to the point of discontinuity, but the value of the 
overshoot amplitude is independent of N . In approximating a perfect low pass filter transfer 
function, the deviations from the theoretical values near the cutoff frequency are usually much 
larger than can be allowed. 


To avoid the sharp cutoff overshoot, instead of making the function zero for all values of 
P > p , it can be continued by a sine function which has the same value and the same derivative at 
P = p as the transfer function and, together with its derivative, becomes zero for a specified 
value of p. Instead of using p directly, however, it is more convenient to use a parameter h p , 
of magnitude corresponding to the change in P during 1/4 cycle of the sine termination function. 
If the ratio of the change in p to h p is included in the argument of the sine function, it forces 
both the termination function and its derivative to have the necessary values at their end points. 
The geometry of the sine termination is shown in Figure 2. Throughout the remainder of this 
discussion we shall simplify the notation for the parameter to h, but it should always be taken to 
mean the h obtained when p is used for the normalized frequency. 


In terms of h, a function which will produce the desired termination effect is 

A(p) = y Q 

/, . 77 Po -P\ 

= y 0 


(26) 
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where y 0 is the amplitude of the sine function 
and p 0 is the value of p at the point where the 
function has the value y 0 . The relationship of 
y 0 and p 0 to h is also shown in Figure 2. As 
may be seen, the quantities y 0 and h will de- 
termine the geometry of the cutoff. The param- 
eter h will permit variation of the slope of the 
sine termination. 

To design a filter with a sine termination, 
h must be as small as possible but such that 
the actual frequency response of the filter does 
not depart from the theoretical response by 
more than a permissible tolerance. (As h ap- 
proaches zero, the filter approaches a sharp 
cutoff filter.) In Figure 3 we see the improve- 
ment offered by sine- terminated filters over 
sharp cutoff filters designed for the same cut- 
off frequencies and with the same number of 
weights. 

So, instead of letting the transfer function 
go to zero immediately for all p > p, we set it 
it equal to A(p) in the range P < p < p 0 + h . 
Then we can write 

f p 

W k = I (1) cos 7rkpdp 

Jo 



FREQUENCY NORMALIZATION PARAMETER, p = f/f c 

Figure 2— Geometry of the sine termination function 
A(p) which is used to provide smoother cutoff for low 
pass filter frequency responses. 


f o+h / . 7 tP-po\ r 1 

y 0 II -sin— 1 cos 77kpdp + (0) cos 77 -kpdp 

V h ' Jp 0 +h 


= w k (0) + w k (a) , 


(27) 


where the first integral is the weight computed for sharp cutoff, and the second is the sine termi- 
nation correction. 

We have already seen that if k = 0 , W k (0 > = W 0 (0 > = p. Similarly, since 


K a) 


= y 0 


p p 0 +h / 

. . ttP-Po\ 

l ( 

Sln 2 h ) 


COS 77kpdp, 
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RELATIVE RESPONSE, W(p) 



FREQUENCY NORMALIZATION PARAMETER, p = f/f c FREQUENCY NORMALIZATION PARAMETER, p - f/f c 


Figure 3~Marked contrast between sharp cutoff and sine-terminated approximations to an ideal filter with low cut- 
off at p = 0.2 is illustrated by the frequency response of two filters with N = 20 and N = 50, respectively. The 
approximation is improved by use of the larger filter. 

we have for k = 0 that 

, , r p 0+h r Po+h ttP-Po 

w 0 (a) = y 0 dp - y 0 sin -2 — — dp 

Jp Jp 



(28) 


(29) 


13 


Atp = p we want A(p) = T(p) and A'(p) = T'(p), so we choose A(P) = 1.0 and A'(P)= o. 

A'(P) = 0 requires that p 0 - P = h , or that p 0 = P + h . By substitution the expression for A(P) = l 
then yields y 0 = 1/2. Using these values for y 0 and p 0 we have finally from Equation 28 that 

W 0 (a> = h , 


and hence that 


w 0 = w 0 <°> +w 0 < a > = p + h. (30) 

In a similar way, if one works on through from Equation 27 (see Reference 8), one obtains for 
k / o that 


m _ fcos 77khl Tsin 77k (P + h)l _ _ n , . Sin 77k (P + h) 

Wk = L ^ J ( ’ h) ^ • (3D 

Tables of F(k,h) may be computed independently of any individual filter and then they will be 
available for particular applications. It will be found from the expression for F(k,h) in Equa- 
tion 31 that an indeterminate form is obtained whenever kh = 1/2. I/Hospitars theorem may be 
used to evaluate the expression in that situation, and it is found that F(k,h) = 0.78540 for kh = 1/2. 

One further correction may be added to the weights in order to normalize the gain to 1.0 at 
p = 0 . Let the value of the k th weight obtained from Equation 31 be designated by L k . Then 


A = 1 - 



and the corrected weight is given by 


Wk=L) < + 2 nTT- (32) 

Once the weights have been computed using Equations 30, 31, and 32, the gain or frequency response 
is easily computed using Equation 10: 


W(p) = W 0 + 2 


E 

k=l 


W k COS 77kp . 


We now have all the formulas necessary for the design of least squares- approximated low 
pass filters. In each case the parameter h can be chosen such as to tailor the cutoff of the filter 
to the specific needs involved. We have already given some indication of the fact that h is sen- 
sitive to the number of weights in the filter, and progressively larger values of h are necessary 
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as one goes to progressively smaller filters (h a 1/N) in order for the sine termination to be 
effectively and efficiently accomplished. This would suggest that a certain amount of experimen- 
tation would be necessary to reveal the optimum value of h for a particular filter. This is the 
value below which the values of the function near the cutoff depart from the theoretical values by 
an unacceptable margin due to the size of the overshoot, and above which the termination gets con- 
tinuously smoother but the accompanying increase in the half-width of the main lobe brings down 
the precision. The effect of filter size n on the width of the main lobe is shown in Figure 4. 

From a study of the cutoff characteristics of a number of filters for various values of h, it 
has been found that at h ~ l/N (if one is using the p notation) the terminal oscillations have almost 
completely disappeared. Although a slightly greater degree of smoothness is obtained as one goes 
to still larger values of h, the predominant effect is merely the broadening of the pass band. On 
the other hand, as one employs progressively smaller values of h, the terminal oscillations grow 
rapidly more significant until the large excursions characteristic of sharp cutoff are obtained. 

Thus to select the proper filter one must compromise between bandwidth and termination 
smoothness, and the limiting factors in the compromise are the minimum cutoff slope that will be 
acceptable for the task at hand and the smallest main-lobe to side-lobe ratio that can be tolerated. 
Once the computation of filter weights and corresponding frequency responses has been pro- 
grammed, it is a simple matter to generate a family of response curves for a particular n and 
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FREQUENCY NORMALIZATION PARAMETER, p = f/f c 

Figure 4— Ultra low pass (cutoff ratio p = 0) filter main lobes for selected values of N. The increase 
in sharpness of a filter due to increasing its size is illustrated. 



various values of h to facilitate the final filter selection. Figure 5 illustrates the three cases of 
h = l/N, 1/2N, and 2/N for a filter with p = 0 and N = 100. Figure 6 is a graphical representation 
of the optimization problem. 

On the question of best filter size, we see that a larger value of n permits sharper sine 
termination, and the undesired frequencies are eliminated more efficiently. Another advantage is 
that the effect of an erroneous input point is much less for a larger filter. However, larger filters 
require longer unbroken runs of data. 

Earlier in this discussion (Equation 1) we described how the filtering is accomplished by 
sliding the filter along the data, applying it to 2N + 1 data points to produce the filtered equivalent 
of the data point lying at the center of the filter (at w 0 ) and then moving each weight to the next 
point in the series and repeating the process. It can be seen that this requires that one have at 
least N data values both preceding and following the time range of interest in order to get the 
filtered equivalents of all data points in that range. Thus a limited number of prior points may 
dictate that a smaller filter be used, or else a scheme may be employed (Reference 8) which in- 
volves the use of progressively larger filters as one gets further into the run of data and con- 
versely near the end. In either event, one steps down to less precisely filtered data throughout 
all or at least part of the range of interest, and this may not be acceptable. In that case one must 
go ahead and use the larger filter and settle for fewer output values. 

It may so happen that the characteristics of the filtering job to be done will suggest particu- 
lar filter size as being most convenient. If it then turns out that this filter will be precise enough 



Figure 5-Dependence of low pass filter cutoff charac- 
teristics on the sine-termination parameter h, where 
N = 100. 
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Figure 6-Regions of significance in filter response opti- 
mization for a given value of N. This is a general illus- 
tration of the design problem created by the morphologi- 
cal changes that occur in the filter frequency response 
when h is varied as in Figure 5. The majority of filter- 
ing problems are perhaps best accommodated by a filter 
with h slightly less than l/N where the oscillations are 
still not too large and the cutoff is reasonably sharp. 
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for the task or will give the desired effect and that there are sufficient data values at each 
end of the run to be used, then there is no selection problem. 

For one application of a low pass filter it was required to produce one which could be used to 
compute hourly averages with the effects of high frequency fluctuations removed. The input con- 
sisted of 2.5 minute surface magnetic field values (H component). There are 25 such values 
inclusive to each hour (0-24), so the natural choice for the task was a 25-point filter (N - 12). 
One-hour periodicities have a corresponding value of p = 0.083 . If this were used as the cutoff 
value, then the slow cutoff of this relatively small filter would allow a considerable fraction of 
higher frequency components to get through. This problem was circumvented by choosing the cut- 
off ratio to be p = o . In this way the slow cutoff itself formed a low pass band which had a gain 
of 0.5 at p = 0.08 (the value of h used), and zero at p = 2h = 0.16 . A filter for which P = o is 
called an ultra low pass filter and essentially gives the trend of the input function. The response 
of the filter used is given in Figure 7, and the weights are listed in Table 1. 



Figure 7— Frequency response of the ultra low pass filter 
designed to produce weighted hourly averages from 2.5 
minute data (i.e., the weights were used in the aver- 
aging of the data for each hour to remove any effect on 
the average due to high frequency components). 


Application of a filter by sliding it along 
the input data gives a "running average" of that 
data. To obtain hourly averages of 2.5 minute 
data, the filter was applied in turn to the 25 
data points inclusive in each hour to produce 
the filtered equivalent of the center point value 
or average for that hour, and then the entire 


Table 1 

Low Pass Filter Weights* for 25-Point Filter 
(N = 12) with h = 0.08. 


W 0 

= 

0.07949 

w t 

= 

0.07817 

W 2 

= 

0.07434 

W 3 

= 

0.06828 

w 4 

- 

0.06046 

W s 

= 

0.05146 

W 6 

= 

0.04189 

w 7 

= 

0.03239 

W 8 

= 

0.02350 

W 9 

= 

0.01566 

W 10 

= 

0.00919 

Wn 

= 

0.00421 

w 
vv 12 

= 

0.00071 


* Because all filters for which weights are given in 
Tables 1-5 are symmetrical, only the weights corre- 
sponding to positive indices are tabulated. Thus in 
applying the full filters W(“k) = W(k) is used. 
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TIME (hours) 

Figure 8— Smoothing of magnetic data provided by the 25-point 
ultra low pass filter (Figure 7) by running it along the data. 
The 2.5 minute output data represents the trend of the input 
data. 


filter was moved over one hour and ap- 
plied again. The values obtained in this 
way are the same as every 25th point on 
the curve described by the 2.5 minute 
values obtained by sliding the filter along 
the run of data point by point. An example 
of the effectiveness of this filter for 
smoothing out high frequency fluctuations 
when applied by sliding it along a run of 
2.5 minute data is given in Figure 8. 


Ultra low pass filters may be labeled according to cutoff point, sine termination char- 
acteristics, and size by means of the following scheme. A filter which is classified as a paabbcc 
filter is one with cutoff aa = 100P, sine termination parameter bb = 100 h p , and size cc = N. 
The use of P in front instead of p will indicate that the values for cutoff and h are given in 
terms of the frequency ratio r = f/f s = p/2 . 


THE CHEBYSHEV APPROXIMATION TO T(f) 

FOR A LOW PASS FILTER 

Another mathematical approach to the computation of filter weights is through the use of 
Chebyshev polynomials (Reference 9). The Chebyshev polynomials are defined by 

T (x) = cos 

m v 


= cos (m arc cos x) . 


The Fourier expressions for orthogonality are 


f 


cos m 0 cos n 8 d 9 


{ 0 f or m / n , 
7 t/ 2 for m = n, 

7 t for m = n = 0, 


and 


N-l 



cos m6 . cos n& . 
j j 


{ 0 f or m / n . 
N/2 for m = n . 

N f or m = n ~ 0 . 


(33) 


(34) 
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In terms of the Chebyshev polynomials these may be written 


i 


T (x) T (x) 


dx 




C 0 for m / n , 

< tt/ 2 f° r m ~ n , 
j for m = n - 0 , 


(35) 


and 


N-l 


i = 0 


r n (x 3 ) 


C 0 for (m / n) , 

< N/2 for (m = n) , 

[ N for (m - n - 0) . 


That the T m (x) are polynomials is shown by the following: De Moivre's theorem states that 


cos m 0 + i sin m 6 = (cos 6 + i sin d) m . 


Expanding the right hand side, taking the real part of both sides, and replacing even powers of 
sin 6 by 


(sin 2 0) k = (1-cos 2 0) k , 


it may be established that cos nd is a polynomial of degree n in cos 6 . However, since 
cos(arc cos x) = x , then T m (x) = cos(m cos x) is a polynomial of degree m in x. 

A general property of orthogonal polynomials is that they are easy to compute and to convert 
to a power series form. From the definition for the Chebyshev polynomials one obtains 


T 0 (x) = l (l). 

T i (x) = X (cos d ) , (36) 

T 2 (x) = 2x 2 - 1 (cos 20) , 

and so on. Another property is that orthogonal polynomials satisfv a three-term recurrence re- 
lation. In the case of the Chebyshev polynomials this relation is 

T n+1 (x) = 2xT n (x)-T n _ 1 (x). (37) 

Finally, we have that |T m (x)| < 1 for |x| < 1 . 

According to the minimax principle, Chebyshev approximations are associated with those 
approximations which minimize the maximum error. Least squares approximation keeps the 
average square error down, but in so doing isolated extreme errors are allowed. Chebyshev ap- 
proximation keeps the extreme errors down but allows a larger average square error. 
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For a symmetrical filter, 


W (f) - W 0 + 2 y Wj cos 277 f j At . 
j=i 

If we let x = cos 6 = cos 7rf At = cos [( 77 / 2 ) f/f c ] = cos (77p/2) , then W(p) is a polynomial of degree 
2m in x . Now also let x = z/z 0 . We have then that T 2m (z) is a polynomial of degree 2m in 
z = xz 0 . By extending the Chebyshev range beyond one then 

|T 2m (z)| > 1 for |z| > 1. 

Because of the similarity of analytical form we equate 

W(p) = T 2m (xz 0 ). (38) 

Now to approximate the frequency response with T 2m (z), it is necessary to adjust T 2m (z) to give 
W(p) = 1.0 at p = 0. For f/f c = p = 0, x = cos(77p/2) = 1 . But x = z/z 0 y so that x = 1 requires 
z = z 0 . Hence T 2m (z 0 ) corresponds to p = 0 . Let T 2m (z 0 ) = r. Then W(p) = T 2m (z)/r will give 
the normalized frequency response (Figure 9). At z = 1 , 

T 2 m ( Z ) =T 2m( 1 ) “ COS f 2m 3rC COS C 1 )^ 

= cos 4m77 
= + 1 . 


Hence, as may be seen in Figure 9, r is the ratio of the main lobe to the amplitude of the side 
lobes. The maximum precision of a particular filter of size N will be realized if the weights W k 

are determined in such a way that the ratio r 
is maximized and the width of the main lobe is 
minimized. 


Now we must investigate the zeros of 
f 2 m ( z ) * Writing 



T 2 m ( XZ o) = COS 2mn » 


Figure 9— Filter frequency response approximation by the 
Chebyshev geometry when the Chebyshev range is ex- 
tended beyond unity (T m (z) > 1 for Z > 1). 


where n = arc cos xz 0 , it is seen that for 
T 2m ( xz 0 ) ~ 0, we must have 


2mn. = k7 t + — . 

k - 2 


By writing 
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2mn, 

k 


(2k - 1) 
2 


then 


2m arc cos x fc z Q = (2k- 1) — , 


or 


arc cos x k z 0 = 


(2k -1) 
4m 


77. 


But 


7T 77 

cos 2 F = cos I p k * 


Therefore 


arc cos 


( z 0 cos 7 Pk ) 




or 


2 o cos ^ P k = cos 


(2k -1), 


( 39 ) 


This will give the values of p k at which the various zeros occur. Most importantly, the first 
zero, which falls at the half-width Pj of the main lobe, is given by 


z o cos ^ Pi = cos (40) 

Hence for any filter size m = N and desired main lobe halfwidth p,, one can compute the corre- 
sponding z 0 from 


4N 


77 

cos T p i 


( 41 ) 


By knowing z Q one can then compute r = T 2m (z 0 ). Finally one can make use of the relation 


W(p) - 


^•2N ( Z ) 


r . 


to compute both the filter weights and the associated frequency response profile. 
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As a brief example, we shall perform the calculations for the weights of a 7-point low pass 
filter (N = 3). Now 


W (p) = w 0 + 2 W. cos 2j 6, where 6 - \ p 

j=i 


= W 0 + 2Wj cos 26 + 2W 2 cos 4 £ + 2W 3 cos 6£. 
Let x = cos 6 . From the recurrence relation we obtain 


(42) 


cos 26 - T 2 (x) = 2x 2 - 1 , 


and 


cos 4 6 - T 4 (x) = 8x 4 - 8x 2 + 1 , 


cos 6 6 = T 6 (x) = 32x 6 -48x 4 + 18x 2 - 1 . 


(43) 


Substituting from the relations in Equations 43 into Equation 42 gives 


W (p) = 32I 3 X 6 + (8 1 2 -48I 3 )x 4 + (2Ij -8I 2 + 18I 3 )x 2 + (I 0 -I 1 + I 2 -I 3 ), 


(44) 


where 


I 3 = 2W 3 , I 2 = 2W 2 , I x = 2W X , and I 0 = W 0 . 


But, as we have seen. 


T 2N ( xz o) _ T 6 ( xz o> _ 32 z o x 6 - 48 z 4 X 4 + 18 z 2 X 2 - 1 ^ 

r r r 

We can now equate T 6 (xz 0 )/r = W(p) by powers of x. This gives 

w 3 = Zg/2r, 

w 2 = 12W 3 - 6zJ/2r, 

Wj = 4W 2 - 9W 3 + 9zg/2r, 

and 

W 0 = 2(Wj - W 2 + W 3 ) - 1/r. (46) 

If, for example, we choose p t =0.5, then for N = 3 we get from Equation 41 that z 0 = 1.37. 
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Then 


r =T 6 (1.37) = 92, 


and from Equation 46 it is found that 
W Q = 0.2509, 

Wj = 0.2113, 

W 2 = 0.1218, 


and 



FREQUENCY NORMALIZATION PARAMETER, p = f/f c 


W 3 = 0.0415 . 

Figure 10— Frequency response of a crude Chebyshev 
low pass filler (a) and the effect produced by applying 
the shifting theorem to that same filter (b). 

The frequency response T 6 (xz 0 ) of this filter 
is shown in Figure 10(a). As can be seen, this 

filter cuts off too slowly to be useful except for some types of smoothing, but it illustrates the 
principle of computing filter weights using Chebyshev polynomials. 


SCALING FILTERS AND SHIFTING FOR BANDPASS RESPONSE 

Whereas the computation of the weights for a large filter using the least squares approxima- 
tion method is no more of a problem than for a small filter if one has automatic computational 
facilities available, the computation of a large Chebyshev filter is extremely complex due to the 
size of the polynomial involved. One way around this is to take a smaller filter and apply scaling 
and interpolation to produce a larger filter. 

Filter scaling is accomplished in the following manner. Suppose we have 


W 2 (t) = 


m 

L 

k =— m 


W k S(t + kAt). 


(47) 


Scaling by a factor a has the effect 


W 2 (t) = 


m 

L 

k-_m 


\ S (a t + kAt) = Wj (at). 


(48) 
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Nonzero weights exist for t = ±kAt/a due to the delta function. Correspondingly, when we 
Fourier transform these functions we obtain 

m t 

Y^(w) = w k S(at + kAt)e- i27rft dt 

J - 00 


or 


= 2 W k cos 


( 49 ) 



COS 



i sin 



- Yj (f/a). 


(50) 


Since nonzero weights correspond to 


_t_ 

At 



(51) 


scaling a filter by a > l means a contraction in the time domain and a corresponding expansion 
in the frequency domain. Conversely, use of a < 1 results in an expansion in the time domain 
and a corresponding contraction in the frequency domain. 


Thus application of filter weights to every m th input point has the same effect on the output 
data as if the sampling frequency had been divided by m . Likewise one can scale a smaller filter 
so as to effect a contraction in the frequency domain until the desired cutoff point is reached, and 
then can use interpolation to increase the number of weights until the desired filter size has been 
reached. Figure 11 shows the response of a 201-point filter (N = 100) computed in this manner. 
In comparison, the figure also shows the response of a least squares filter of the same size. The 
weights are tabulated in Tables 2 and 3. 


Much of the power of the numerical filtering technique comes from the possibility of being 
able not only to low pass filter data, but also to attenuate all frequencies outside a particular fre- 
quency band of interest and hence bandpass filter data as well. The simplest way to accomplish 
this is to shift a low pass filter from f = o to f = q, the central frequency of the band of interest. 
The sharpness of response of the low pass filter determines the effective width of the bandpass 
window. 
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Wj (f ) = f W(t)e" i27Tft dt, 

J -CO 

a shift by r in time gives 


(52) 


W_ (f) — W(t +T)e' i2,rf <* +T > d(t +r) 


£ 

r 


W(t)e _i2,rft e iivtT dt 


e -27rfT W(t)e~ i2,Tft dt 

~ —CO 

e -i 2rrfr ffj (f) . (53) 



Figure 11— Example of response of (a) least squares and (b) 
Chebyshev approximations of ultra low pass filters (201- 
pomt, N = 100). The least squares filter was derived with 
h = 0.01, i.e., h = l/N. If the value of h used in the 
least squares approximation had been slightly less than l/N, 
the response would have been nearly that of the Chebyshev 
filter. 


Hence a shift by r 


in time leads to a shift by e i27rfT 


in the frequency domain. Since 


and 


W(f) = J W(t)e- i2wft dt 


W(t) = 


f W(f)e +i27rft df 
J —00 


are identical except for a change in sign, if one desires to shift by q in frequency, one must multi- 
ply by e i27Tqt in time. In order to achieve a bandpass filter at f = q that has a symmetrical 
transfer function, negative frequencies must be considered as well and hence a low pass filter 
must be shifted by ±q. 

We define 


W BP (f) = 8 [W LP ( f + q) + W LP (f -q)]- (54) 

By the shifting theorem described above this gives 

w 813 (t) = (e i27rqt + e‘ i27rqt ) W LP (t) 

= 2 cos 277qt W LP (t). (55) 
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Table 2 


Low Pass Filter Weights 

for 201 -Point Expanded Chebyshev Numerical Filter. 


W ( 

0) 

= 

0.0085210 

W< 

51 ) 

= 

C. 0050260 

W( 

1) 

= 

0.U085200 

w< 

52) 

= 

0.0049170 

W( 

2) 

= 

0.0085150 

W( 

53) 

= 

C. 0048070 

W( 

3) 

= 

0 • 0C 65070 

Vt ( 

54) 

= 

0.0046980 

W( 

4) 

= 

0.0084950 

w< 

55) 

= 

C. 0045880 

W( 

5) 

= 

0.0084800 

W( 

56) 


0.0044790 

W( 

6) 

= 

0.0084620 

W( 

57) 

= 

0.0043690 

W( 

7) 

= 

0.0084410 

j W{ 

58) 

= 

0.0042600 

W( 

8) 

= 

0.0084160 

W{ 

59) 

= 

C. 0041520 

W( 

9) 

= 

0.0083880 

W( 

60) 

= 

0.0040440 

W( 

10) 

= 

0.0083580 

w ( 

61) 

= 

0.0039360 

W( 

11) 

= 

0. 0083240 

W( 

62) 

= 

0.0038290 

w( 

12) 

= 

0.0082860 

K ( 

63) 

= 

C. 0037220 

W( 

13) 

= 

0.0082460 

W( 

64 ) 

= 

0.0036170 

W( 

14) 

= 

0.0082030 

W( 

65 ) 

= 

C. 0035110 

W( 

15) 

= 

0.0081570 

W( 

66) 

= 

0.0034070 

W( 

16) 

= 

0.0081070 

W( 

67) 

= 

0.0033040 

W( 

17) 

= 

0.0080550 

W( 

68 ) 

= 

0.0032020 

W ( 

18) 

= 

0 . 008 C 000 

w ( 

69) 

= 

0.0031010 

W ( 

19) 

= 

0.0079420 

W( 

70 

= 

0.0030000 

W ( 

20) 

= 

0.C07882C 

W( 

71 ) 

= 

C. 0029C2 0 

W{ 

21) 

= 

0.0078190 

W{ 

72) 

= 

0.0028040 

W( 

22) 

= 

0. C077530 

W{ 

73 ) 

= 

C. 0027070 

W( 

23) 

= 

0.0076840 

W( 

74) 

= 

0.0026120 

w ( 

24) 


0.C076130 

W( 

75) 

= 

C. 0025190 

w< 

25) 

= 

0 • C075400 

K( 

76) 

= 

0.0024260 

W( 

26) 

= 

0.0074640 

w ( 

77) 


0.0023350 

W( 

27) 

= 

0.0073860 

k ( 

78) 

= 

0.0022460 

W( 

28) 

= 

0.C073060 

W{ 

79) 

= 

C. 0021580 

W ( 

29) 

= 

0.0072230 

K{ 

80) 

= 

0.0020720 

W( 

30) 

= 

0.C071390 

W( 

81 ) 

= 

0.0019870 

W( 

31) 

= 

0. CO 70520 

W ( 

82) 

= 

0.0019040 

W( 

32) 

= 

0.0069640 

W( 

83) 

= 

C. 0018230 

W( 

33) 

= 

0.0068730 

W( 

84) 

= 

0.0017430 


34) 

= 

0.C067810 

W( 

85) 

= 

0.0016650 

W( 

35) 

= 

0.0066870 

W( 

86) 

= 

C. 0015890 

W( 

36) 

= 

0.0065920 

W( 

87) 

= 

0.001515C 

W ( 

37) 

= 

0.0064950 

W( 

88 ) 

= 

0.0014420 

W( 

38) 

= 

0.0063960 

W( 

89) 

= 

0.0013710 

W( 

39) 

= 

0.0062970 

W( 

90) 

= 

0.0013020 

W( 

40 ) 

= 

0.0061960 

M 

91) 

= 

C. 0012350 

w< 

41) 

= 

0.0060940 

W( 

92 ) 

= 

0.0011700 

W( 

42) 

= 

0.0059900 

W ( 

93 ) 

= 

0.0010170 

W( 

43) 


0.0058860 

41 

94) 

= 

0.0009260 

W( 

44) 

= 

0.0057810 

VI ( 

95) 

= 

0.0008960 

W{ 

45) 

= 

0.0056750 j 

V ( 

96) 

= 

0.0009270 

W( 

46) 

ss 

0.0055680 

U( 

97) 

= 

0.0010190 

W( 

47) 

= 

0.0054610 


98) 

= 

0.0011730 

W( 

48) 

= 

0.0053530 

w< 

99) 

= 

0.0013870 

W ( 

49) 

= 

0.0052440 

W( 100 ) 

= 

0.0016630 

W( 

50) 

= 

C. C051350 







Table 3 


Low Pass Filter Weights 

for 201- Point Least Squares Numerical Filter (P = 0, h = .01) 


to ( 

C ) 

= 

C.CC99368 

tot 

51) 

= 

0.CC47668 

M 

1 ) 

- 

0.CC99342 

to t 

52) 

= 

G.CC46377 

tot 

2 ) 

= 

O.OC99265 

tot 

53 ) 

= 

C.CC44893 

tot 

3 ) 

= 

C.CC99136 

tot 

54 ) 

= 

C.CC43416 

to ( 

4 ) 

= 

0.CC98956 

tot 

55) 

= 

C.CC41949 

to ( 

5 ) 

x 

0.CC98725 

to ( 

56) 

= 

C.CC4C493 

tot 

6) 

x 

G.C098443 

tot 

57) 

= 

C.CC39C49 

tot 

7) 

st 

0.CC9SU1 

to ( 

50) 

X 

0.CC37619 

v% c 

8) 

- 

0.CC97728 

to ( 

59 ) 

X 

C.CC362C3 

to ( 

9 ) 

- 

0.CC97296 

tot 

60 

= 

C.CC348C3 

ini ( 

1C ) 

X 

0.CC96815 

tot 

6 1 ) 

X 

C.CC3342C 

to ( 

n ) 

X 

C.CC96286 

tot 

62) 

= 

0.CC32C55 

Vs ( 

12 ) 

x 

0.CC957C6 

tot 

63) 

X 

C.CC3C7C9 

tot 

13) 

= 

C.CC95C84 

tot 

64) 

X 

C.CC29382 

Vs ( 

14 ) 

x 

0.CC944 13 

vs t 

65) 

X 

0.CC28077 

Vs ( 

15 ) 

x 

O.CC 93697 

tot 

66) 

= 

0.CC26793 

Vs ( 

16) 

= 

0.CC92936 

tot 

67) 

X 

C.CC25532 

Vs ( 

1 7 ) 

= 

0.CC92132 

to t 

68) 

= 

C.CC24255 

Is ( 

lb ) 

X 

C.CC91285 

tot 

69) 

= 

0.CC23061 

V* ( 

IS ) 

= 

0.CC9C396 

to t 

7 C ) 

X 

C.CC21693 

Vs l 

2C ) 

= 

O.CCd94fc7 

tot 

7 1 ) 

= 

C.CC2C73C 

Vs ( 

21 ) 

= 

0.C066498 

to ( 

72) 

= 

O.CC 19593 

Vs ( 

22 ) 

= 

O.C08 74 9 C 

tot 

73) 

= 

C.CC18483 

Vs ( 

23 ) 

= 

0.CCS6446 

tot 

74 ) 

= 

C.CC174CC 

tot 

24 ) 

= 

0 • C 0 6 5 36 6 

to ( 

75) 

X 

C.CC16345 

Vs { 

25) 

= 

C.C064251 

tot 

76) 

X 

C.CC15318 

Vs ( 

26 ) 

X 

O.CCd3lC3 

to ( 

77) 

X 

C.CC14319 

Vs ( 

2 7 > 

- 

O.CCd L923 

tot 

76) 

X 

C.CC13349 

Vs ( 

2b) 

X 

O.CC8C712 

to { 

79 ) 

X 

C.CC124C8 

Vs < 

2 9 ) 

X 

O.CC7947? 

to ( 

8 0 ) 

X 

C.CC 11497 

Vs ( 

3 6 ) 

= 

0.CC782C4 

to 1 

e i ) 

X 

C.CC1C615 

Vs { 

3 i ) 

X 

G * C C / 6 9 l i 

w( 

62 ) 

X 

C.CCC9762 

Vs { 

32 ) 

= 

C.CC 7 5 592 

tot 

83) 

= 

C.CCC894C 

tot 

33 ) 

X 

O.CC 7 4 2 5 C 

tot 

64 ) 

X 

C.CCC8147 

Vs < 

34 ) 

X 

O.CC 72866 

tot 

85 1 

X 

C.CCC7363 

Vs ( 

3 6 ) 

= 

O.CC 715C2 

to ( 

86 ) 

= 

C.CCC6649 

vs { 

36 ) 

X 

O.CC 7C ICC 

to t 

6 7) 

X 

0.CCC5 745 

Vs ( 

3 7 ) 

X 

0.CC6668C 

tot 

6 8 ) 

X 

0.CCC5271 

Vs ( 

3 b ) 

X 

C.CC 6 7244 

to t 

69) 

X 

C.CCC4625 

Vs { 

3 S ) 

X 

G.CC65795 

tot 

90 ) 

= 

C.CCC4CC9 

tot 

4C ) 

X 

0.CC64333 

to { 

5 l ) 

X 

C.CCC3421 

Vs l 

4 1 ) 

X 

C.CC6286C 

to « 

92) 

= 

C.CCC2862 

to ( 

4 2 ) 

X 

0 . C C fc 1 3 7 8 

tot 

9 3 ) 

X 

C.CCC2331 

to ( 

4 3 ) 

X 

o.cc59eee 

tot 

94 ) 

X 

C.CCC 1828 

to ( 

44 ) 

X 

O.CC 6 6392 

to ( 

95 ) 

X 

C.CCC1352 

to ( 

4b ) 

X 

0.CC56891 

to ( 

96) 

X 

C.CCCC9C3 

to ( 

4 7 ) 

X 

C.CC5536 / 

vs { 

97) 

X 

C.CCCC461 

to ( 

4 7 ) 

X 

C.CC53681 

tot 

98 ) 

X 

C.CCCCC85 

to { 

4 h ) 

X 

C.CC 52376 

tot 

9 9 ) 

*• 

-C.CCCC2R6 

tot 

4 S ) 

X 

C.CC5C972 

to { 

ICO 

= ■ 

-C.CCCC6 32 

to ( 

51 ) 

- 

C.CC 4 9 366 







Thus to compute the k th weight of a bandpass filter centered at f = q , one must operate in the 
following way on the k th weight of the low pass filter to be shifted: 

W® p = 2 cos 77 kq/f c W]/\ (50) 

The response of a bandpass filter obtained by performing this operation on each weight of the 
N = 3 low pass Chebyshev filter computed previously as an example is shown by curve B in 
Figure 10. 

In order to filter out the diurnal component of geomagnetic data, bandpass filters were con- 
structed by shifting the expanded Chebyshev filter (N = 100) shown in Figure 11. The first five 
harmonics of the total diurnal component have periods of 24, 12, 8, 6 and 4.8 hours. Thus the LP 
filter was successively operated upon to give five BP filters which peaked at p = 1/12, 2/12, 3/12, 
4/12, and 5/12. Power spectral analysis performed on a number of typical runs of H-component 

data during the IGY (Ness, 1962*) has revealed 
that one need not be concerned with more than 
the first five harmonics of the diurnal com- 
ponent. This may also be ascertained by in- 
spection of the relative amplitudes of the 
harmonics themselves. Figure 12 shows the 
first five harmonics at a middle latitude ob- 
servatory for a period during which the field 
was disturbed and thus contained harmonics of 
enhanced amplitude. 

The individual harmonics shown in Fig- 
ure 12 were isolated from geomagnetic H com- 
ponent data by using each of the five BP filters 
individually on the raw data. Although this 
sort of analysis reveals the relative amounts 
of energy present in the different harmonics 
and shows modulation due to time variations in field strength, our primary intent was to develop a 
filtering process that would isolate the total diurnal component so that it could be subsequently 
subtracted from the raw surface data. The most efficient way to achieve this effect was to combine 
the five BP filters in such a way as to produce one resultant BP filter. 

Each filter, when applied to the raw data, produces as its output the particular harmonic com- 
ponent which it was designed to pass. In our case 
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Figure 12-The first five harmonics of the diurnal com- 
ponent of the horizontal magnetic field at Fredericksburg, 
Virginia during the week of geomagnetic disturbances of 
December 2, 1963. 


*N. F. Ness, private communication. 



0,(t) =F„(t) 
° 2 (t) = F ia <t) 
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( 57 ) 


To a good approximation the diurnal component is a linear function of only these five harmonics, 
i.e., 


D = F 24 (t) +F I2 (t) +. . . *F 4-8 (t). 


(58) 


Since this is equivalent to 


O f (t) = 0,(t) + o 2 (t) + . . . + o s (t). 


(59) 


and since we have for an input function i(t) that 


0,(t) 


■L 


I(t-j 


and similarly for the other components, then we can write 


(60) 


o f (t) 


i ( t - ) At > wj 1 


L 


I(t - j At) Wj 2 * * + . . . 


+ 


Z>- 


j At) Wj 5 


I (t - jAt) (Wj + W2 + . . . 4- W j) 


= D. (61) 

Hence the output of the filter obtained by linearly combining the corresponding weights from each 

separate BP filter is the total diurnal component that we wish to isolate. In general, for the k th 
weight of a BP filter which will pass only a data component which is the sum of m harmonics, we 

simply add: 
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w k = K + w k + • • • + K - 


( 62 ) 


where is the k th weight of the i th harmonic filter. 

The frequency response of the filter which was constructed for isolating the total diurnal 
component is shown in Figure 13, and the weights are listed in Table 4. The data values com- 
prising the output of this filter were subtracted from raw data values of corresponding times, 
producing geomagnetic H-component data free of the diurnal modulation. In Figures 14 and 15 
examples of hourly average data and the corresponding filtered data have been plotted. 

One further application of our filters involved plotting on an expanded scale the raw 2.5 minute 
H component data for the time period covered by magnetic storms. We wanted the diurnal com- 
ponent removed from the data, but to prevent aliasing we had to low pass the data to remove the 
high frequency constituents before application of the bandpass filter. The output of the LP filter 
was used as input to the BP filter, and the output of this second filter was the dirunal component. 

Scaling was used in the application of the BP filter to the extent that the weights were applied 
to 2.5 minute data values that were separated in time by one hour. Since the sampling frequency 
of the data was 24/hour, application of the weights to every 24 th data value was equivalent to 
dividing the sampling frequency by 24, giving the one hour sampling frequency required by the 
design of the filter. The residual that remained when the output of the BP filter was subtracted 
from the raw data was the H-component of the geomagnetic field as described by 2.5 minute data 
points. 


Figure 16 is an example of the results of applying this filtering technique to the data recorded 
at six magnetic observatories during the storm of April 1, 1964. With the diurnal components re- 
moved, the remaining time variations can be 
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Figure 13— Frequency response of filter derived by linear 
combination of five 201 -point bandpass filters obtained 
by shifting the expanded Chebyshev filter band shown in 
Figure 11 to values of p corresponding to the frequen- 
cies of each of the first five harmonics of the diurnal 
component. 


more effectively correlated with the field 
variations observed at a distance from the 
Earth by the satellite. 

Another filter application problem re- 
quired the construction of a bandpass filter 
centered on the 24-hour component and with a 
±4-hour bandwidth, i.e., the passband had to 
be centered on p 24 = 0.0833, with cutoffs at 
p 2 s = 0.0714 and p 20 = 0.1000. It was further 
required that N = 24, or perhaps some larger 
multiple of this if N = 24 did not give satis- 
factory results. 

It is difficult to design a filter of this size 
or smaller with such a narrow passband be- 
cause the wavelength of the oscillations in the 
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Table 4 


Bandpass Filter Weights 

for 201-Point Harmonic Diurnal Component Filter. 


W( 

0) 

= 0.0852100 

W{ 

51 ) 

=-0.0171596 

w< 

1 ) 

= 0.0561956 

W< 

52) 

=-0.0098348 

w< 

2) 

=— 0 . 

W( 

53) 

= 0.0014574 

W { 

3) 

=-0.0290446 

W( 

54) 

= C * 0 C COO 04 

W( 

4) 

=-0.0169900 

Wl 

55) 

=-0.0081082 

W( 

5) 

= 0.C025714 

W( 

56) 

=-0.0089580 

W ( 

6) 

= 0.0000002 

w< 

57) 

=-0.0025596 

W{ 

7) 

=-0.0149180 

w( 

58) 

= C.00CC0C2 

W( 

8) 

=-0.0168320 

W( 

59) 

=-0.0046982 

W( 

9) 

=-0.0049138 

W( 

60) 

=-0.0080880 

W( 

10) 

= 0.0000002 

W( 

61 ) 

=-0.0044546 

W( 

11) 

=-0.0094198 

W( 

62) 

= — 0 . 

W( 

12) 

=-0.0165720 

w( 

63) 

=-0.0021802 

W( 

13) 

=-0.0093316 

Vl( 

64) 

=-0.0072338 

W ( 

14) 

=-0.0000002 

M 

65) 

=-0.0062054 

W( 

15) 

=-0.0047780 

W( 

66) 

=— C. 0CCCC04 

W( 

16) 

=-0.0162138 

w( 

67) 

= 0.0010022 

W( 

17) 

=-0.0142360 

w ( 

68) 

=-0.0064034 

W ( 

18) 

=-0.000000 2 

m ( 

69 ) 

=-0.0105876 

H( 

19) 

= 0.0024082 

w ( 

70 ) 

=-0.0000012 

W( 

20 ) 

=-0.01 57638 

w ( 

71 ) 

= 0.0191398 

W( 

21 ) 

=-0.0266958 

W( 

72) 

= 0.0280400 

W ( 

22) 

=-0.0000008 

w ( 

73) 

= 0.0178558 

W( 

23) 

= 0.0506810 

Mi ( 

74) 

= 0.0000012 

W( 

24) 

= 0.0761300 

w ( 

75) 

=-0.0086004 

W( 

25) 

= 0.0497330 


76) 

=-0.0048524 

W( 

26) 

= 0.0000010 

W( 

77) 

= 0.0007080 

W( 

27) 

=-0.0252172 

w( 

78 ) 

= 0.000C002 

W( 

28) 

=-0.0146126 

W ( 

79) 

=-0.0038138 

W ( 

29) 

= 0.0021900 

w ( 

80) 

=-0.0041444 

W( 

30) 

= 0.0000002 

W( 

81 ) 

=-0.0011642 

W( 

31) 

=-0. 01 2*628 

w { 

82 ) 

=— 0 • 

W ( 

32) 

=-0.013928? 

W( 

83) 

=-0.0020630 

W ( 

33) 

=-0.0040262 

W( 

84) 

=-0.0034860 

W ( 

34) 

=-0.0000002 

M ( 

85) 

=-0.0018844 

W ( 

35) 

=-0.0075670 

W( 

86) 

=— 0 • 

W( 

36) 

=-0.0131840 

w( 

87) 

=-0.0008872 

W ( 

37) 

=-0.0073504 

W( 

88) 

=-0.0028838 

W ( 

38) 

=— 0 * 

W ( 

89) 

=-0.0024232 

W( 

39) 

=-0.0036886 

w ( 

90 ) 

=— 0 • 

W( 

40 ) 

=-0.0123918 

W( 

91 ) 

= 0.0003746 

W ( 

41) 

=-0.0107704 

W( 

92) 

=-0.0023396 

w< 

42) 

=-0.0000004 

W( 

93) 

=-0.0034722 

W( 

43) 

= 0.0017852 

W( 

94) 

=-0.0000006 

W( 

44) 

= -0-. 01 15616 

W( 

95) 

= G. 0059092 

W( 

45) 

=-0.0193756 

W( 

96) 

= 0. 00927C0 

W( 

46) 

=-0.0000014 

W( 

97) 

= 0.0067216 

W( 

47) 

= 0.036018G 

n< 

98) 

= 0. 00CC006 

W( 

48 ) 

= 0.0535300 

w( 

99) 

=-0.0047354 

W( 

49) 

= 0.0345894 

W( 

100) 

=-0.0033264 

W( 

50) 

= 0.0000014 
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Figure 14— Hourly average data representing the varia- 
tions in the H-component of the geomagnetic field at 
seven magnetic observatories from May 9 to June 3, 
1964. The raw 2.5 minute data values from these sta- 
tions were averaged using the weights of the filter shown 
in Figure 7. 



1964 

MO NTH /DAY 5/9 5/14 5/19 5/24 5/29 6/3 

Figure 15— The output data obtained when the data shown 
in Figure 14 are used as inputs to the diurnal component 
BP fi Iter of Figure 13. 


DISTANCE (earth radii) 



Figure 16— Result of applying the diurnal component 
BPfilter to 2.5 minute geomagnetic data from six mag- 
netic observatories. These filtered data show the H- 
component variations that occur at each of the stations 
prior to and during the magnetic storm of April 1, 1964. 


approximating function (the actual frequency 
response) is greater than the width of the 
theoretical band. This prevents a very close 
approximation to the rectangular shape of the 
ideal band. Also in order to approach a step 
function type profile for the low pass filter to 
be shifted, one must employ a cutoff value of 
p > 0. However, when P is increased in an 
attempt to fit the W(p) =1 part of the step, one 
finds that there is an accompanying rapid in- 
crease in the width of the achieved band at the 
base (near the W(p) = 0 level). Some results for 
N = 24, 48 and 72 are shown in Figures 17-19, 
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Figure 17— Frequency response obtained in the attempt 
to approximate a 24± 4 hour ideal passband with a re- 
latively small filter (49 -point, N = 24). 



FREQUENCY NORMALIZATION PARAMETER, p = f/f c 

Figure 18— Improved approximation of the 24 ± 4 hour 
ideal passband obtained by doubling the filter size (97- 
point, N = 48). 


and the corresponding weights are tabulated 
in Table 5. As can be seen, it is only when 
N = 72 is used that a close approximation of the 
desired band is approached. 

CONCLUDING REMARKS 

In summary, an attempt has been made to 
present (1) a brief review of the basic concepts 
of numerical filtering, (2) derivations of formu- 
las which can be utilized to compute filter 
weights and frequency profiles by using either a 
least squares approach or Chebyshev polyno- 
mials to approximate the desired response, and 
(3) several examples of how such filters are 
currently being used in studying transient variations of the earth T s magnetic field. The presen- 
tation has been oriented more toward providing a guide to practical application than a rigorous 
theoretical treatment of the subject. 

For supplemental information on numerical filtering see Fullenwider and McNamee (Refer- 
ence 11), and Martin (References 8 and 12). Arthaber (Reference 13) has considered the related 
topic of spatial filtering. The more general subjects of digital filtering and frequency analysis 
have been discussed in the literature by Salzer (References 14 and 15), Linville and Salzer (Refer- 
ence 16), Hauptsehein (Reference 17), and Langill (Reference 10). 
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Figure 19— Close approximation (compare Figures 17 and 
18) of the 24 ± 4 hour ideal passband obtained by use of 
a 145-point filter (N increased to 72). 
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Table 5 

Bandpass Filter Weights, for 24 ± 4 Hour Passband. 





49 Points 

h = .025 p = .020 

97 Points 

h = .014 p = .016 

145 Points 

h = .012 p = .012 

W ( 0) 


0.08799 

W( 0) 

= 

0.05969 

W126) 

= 

0.00883 

W ( Ui 

= 

0.04769 




W( 1) 

— 

0.08445 

W ( 1) 

= 

0.05738 

W ( 2 7 ) 

= 

0.00492 

Wl 1) 

= 

0.04591 

W l 37 ) 

= 

-0.00423 

W( 2) 


0.07422 

W( 2) 

- 

0.05065 

W ( 28 ) 

= 

0.00 173 

W( 2) 


0.04C72 

W ( 38 ) 

= 

-0.00258 

W ( 3) 

= 

0.05841 

Wl 3) 

= 

0.04011 

W ( 29 ) 


-0.00048 

W( 3) 

= 

0.03255 

W ( 39 ) 

= 

-0.00123 

Wt 4) 


0.03874 

Wl 4) 

= 

0.02670 

WOO) 

= 

-0.00 159 

Wl 4) 

= 

0.02206 

W l 40 ) 


-0.00033 

W ( 5) 

- 

0.01728 

W( 5) 

= 

0.01 160 

W ( 31 ) 

= 

-0.00167 

W( 5) 

= 

0.01013’ 

W { 41 ) 

•= 

0.00006 

Wl 6) 


-0.00380 

W( 6) 

- 

-0.00 18 8 

W ( 3 2 > 

= 

-0.00091 

Wt 6) 

= 

-0.00230 

W ( 42 ) 

= 

-0.00004 

Wl 7) 

— 

-0.02247 

i W{ 7) 

= 

-0.01843 

W ( 33 ) 

= 

0.00037 

Wl 7) 

= 

-0.01424 

W (43 ) 

= 

-0.00056 

W ( 8) 

- 

-0.03713 

' Wl 8) 

= 

-0.03088 

W ( 34 ) 

= 

0.00184 

W( 8) 

= 

-0.02477 

W C 44 ) 

= 

-0.00136 

W( 9) 

- 

-0.04674 

W( 9) 

= 

-0.04029 

W ( 35 ) 

= 

0.00815 

W( 9) 

= 

-0.03310 

W l 4 5 ) 


-0.00225 

WllO) 

- 

-0.05C93 

WHO) 


-0.04604 

W( 36) 

= 

0.00402 

W110) 


-0.03867 

W ( 46 ) 

= 

-0.00307 

W ( 1 1 1 


-0.04599 

Will) 

= 

-0.04 f & 6 

W ( 37 ) 

- 

0.00*21 

W( 11) 

= 

-0.04114 

W ( 47 ) 

= 

-0.00365 

W ( 12 ) 

= 

-0.04477 

W (12 ) 

- 

-0.04587 

W { 38 ) 

- 

0.00365 

W ( 1 2 ) 

= 

-0.04047 

W { 48 ) 

= 

-0.00386 

WI13) 


-0.03655 

W (1 3 ) 

= 

-0.04051 

W ( 39 ) 

- 

0.00237 

W ( 1 3 ) 

= 

-0.03684 

W ( 49 ) 

= 

-0.00363 

W(1M 

~ 

-0.02681 

W (14 ) 

= 

-0.03253 

W ( 40 ) 

= 

0.00^49 

W( 14) 

= 

-0.03072 

W l 50 ) 

= 

-0.00296 

W ( 1 5 ) 

= 

-0.01703 

W ( 15 ) 

= 

-0.02284 

W C 4 1 ) 

= 

-0.00176 

W( 15) 

= 

-0.02271 

W ( 5 1 ) 

= 

-0.00188 

W( 16) 


-0.00846 

W( 16) 

= 

-0.01243 

W { 42 ) 

- 

-0.00412 

W l 1 6 ) 

= 

-0.01 357 

W ( 52 ) 

= 

-0.00049 

W ( 17 ) 

= 

-0.00202 

W H 7 ) 

= 

-0.00 ^30 

W (43 ) 

= 

-0.00631 

W ( 17 ) 

= 

-0.00409 

W 1 53 ) 

= 

0.00106 

W( 18) 

= 

0.00180 

W( 18) 

= 

0.00669 

W (44 ) 

= 

-0.00HO6 

W ( 18 ) 


0.00495 

W ( 54 ) 


0.00262 

W (19 ) 

= 

0.00301 

W { 19 ) 

= 

0.01387 

W ( 45 ) 

= 

-0.00917 

Wl 19) 

= 

0.01285 

W ( 55 ) 

= 

0.00404 

W ( 20 ) 


0.00197 

W ( 20 ) 

= 

0.01882 

W 1 46 ) 

= 

-0.00953 

W ( 20 ) 

= 

0.01910 

W ( 56 ) 

= 

0.00516 

W ( 2 1 ) 

= 

-0.00063 

W ( 2 1 ) 

= 

0.02 139 

W ( 47 ) 

= 

-0.00909 

W ( 2 1 ) 


0.02332 

W ( 57 ) 

= 

0.00589 

W ( 22 ) 


-0.00395 

W ( 22 ) 

= 

0.02169 

W (48 ) 


-0.00792 

W ( 22 ) 

= 

0.02536 

W ( 58 ) 

= 

0.00614 

W ( 23 ) 

= 

-0.00714 

W ( 23 ) 

= 

0.02005 




W ( 23 ) 

= 

0.02527 

W ( 59 ) 

= 

0.00590 

W 1 24 ) 

- 

-0.00947 

W(24) 

= 

0.0169/ 




W ( 24 ) 

= 

0.02328 

W ( 60 ) 

= 

0.00522 




W ( 25 ) 

= 

0.01303 




W l 25 ) 

= 

0.01977 

W (61 ) 


0.00416 










W ( 26 ) 

= 

0.01521 

W ( 62 ) 

= 

0.00284 










W l 27 ) 

= 

0.01011 

W 1 63 ) 

= 

0.00 138 










W { 28 ) 

- 

0.00499 

W ( 64 ) 


-0.00007 










W ( 29 ) 


0.00030 

W ( 65 ) 

= 

-0.00138 










WOO) 

= 

-0.00360 

W (66 ) 

* 

-0.00246 










W { 3 1 ) 

= 

-0.00648 

W ( 67 ) 

= 

-0.00322 










W ( 32 ) 

= 

-0.00823 

W ( 68 ) 

= 
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W ( 33 ) 

= 

-0.00888 

W (69 ) 


-0.00372 

= 

1 









W ( 34 ) 

= 

-0.00856 

W ( 70 ) 

= 

-0.00348 










W 1 35 ) 


-0.00749 

, W ( 7 1 ) 

« 

-0.00300 

1 

i 
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W ( 36 ) 

= 

-0.00596 

W ( 72 ) 


-0.00235 
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conducted so as to contribute ... to the expansion of human knowl- 
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